"""
Natural Convection Demo - Rayleigh-Bénard Convection

This example demonstrates coupled fluid flow and heat transfer with Boussinesq
buoyancy effects. It simulates the classic Rayleigh-Bénard convection problem
where a fluid layer is heated from below and cooled from above, leading to
natural convection patterns.

Key Physics Demonstrated:
- Incompressible Navier-Stokes equations
- Energy equation with convection-diffusion
- Boussinesq approximation for buoyancy
- Temperature-dependent density effects
- Natural convection heat transfer

Problem Setup:
- 2D rectangular domain with heated bottom wall and cooled top wall
- Adiabatic side walls (no heat flux)
- No-slip boundary conditions for velocity
- Rayleigh number Ra = 10^4 to 10^6 range for interesting dynamics
"""

using NSEMSolver
using Printf

"""
    run_natural_convection_demo(; Ra=1e5, Pr=0.71, resolution=32, verbose=true)

Run natural convection simulation with specified parameters.

# Arguments
- `Ra::Float64`: Rayleigh number (default: 1e5)
- `Pr::Float64`: Prandtl number (default: 0.71 for air)
- `resolution::Int`: Grid resolution (default: 32)
- `verbose::Bool`: Enable detailed output (default: true)

# Returns
- `MultiPhysicsResult`: Complete simulation results
"""
function run_natural_convection_demo(; Ra::Float64=1e5, Pr::Float64=0.71, 
                                    resolution::Int=32, verbose::Bool=true)
    
    if verbose
        println("🌡️ Natural Convection Demo - Rayleigh-Bénard Convection")
        println("=========================================================")
        println("Rayleigh number: Ra = $(Ra)")
        println("Prandtl number:  Pr = $(Pr)")
        println("Grid resolution: $(resolution)×$(resolution)")
        println()
    end
    
    # Problem parameters
    L = 1.0              # Characteristic length (m)
    ΔT = 10.0            # Temperature difference (K)
    T_ref = 300.0        # Reference temperature (K)
    g = 9.81             # Gravitational acceleration (m/s²)
    
    # Compute fluid properties from Rayleigh and Prandtl numbers
    # Ra = gβΔTL³/(να), Pr = ν/α
    # This gives us: ν = √(gβΔTL³/(Ra*Pr)), α = ν/Pr
    
    β = 1/T_ref          # Thermal expansion coefficient (1/K, ideal gas approximation)
    ν = sqrt(g * β * ΔT * L^3 / (Ra * Pr))  # Kinematic viscosity
    α = ν / Pr           # Thermal diffusivity
    
    if verbose
        println("Computed fluid properties:")
        println("  Thermal expansion coefficient: β = $(β) 1/K")
        println("  Kinematic viscosity: ν = $(ν) m²/s")
        println("  Thermal diffusivity: α = $(α) m²/s")
        println()
    end
    
    # Base Navier-Stokes solver options
    ns_options = NSOptions(
        N = 5,                    # Polynomial order
        nu = ν,                   # Kinematic viscosity
        cfl = 0.3,                # Conservative CFL for stability
        tfinal = 20.0,            # Final time (multiple convection time scales)
        dim = 2,                  # 2D simulation
        n_block = div(resolution, 5),  # Number of spectral elements
        max_steps = 5000,         # Maximum time steps
        tol = 1e-6,              # Convergence tolerance
        solver = :julia,          # Linear solver backend
        verbose = verbose
    )
    
    # Boussinesq approximation options
    boussinesq_options = BoussinesqOptions(
        reference_temperature = T_ref,           # Reference temperature
        thermal_expansion = β,                   # Thermal expansion coefficient
        gravitational_acceleration = [0.0, -g], # Gravity vector (downward)
        density_reference = 1.0,                # Reference density
        prandtl_number = Pr                     # Prandtl number
    )
    
    # Heat transfer options
    heat_options = HeatTransferOptions(
        thermal_conductivity = 1.0,              # Thermal conductivity
        thermal_diffusivity = α,                 # Thermal diffusivity
        heat_capacity = 1.0,                     # Specific heat
        density = 1.0,                           # Density
        heat_source = nothing,                   # No volumetric heat source
        boundary_conditions = create_rayleigh_benard_bc(T_ref, ΔT),
        temperature_dependent_properties = false,
        convection_scheme = :upwind,             # Upwind for stability
        time_integration = :explicit             # Explicit time integration
    )
    
    # Multi-physics coupling options
    multiphysics_options = MultiPhysicsOptions(
        include_heat_transfer = true,
        include_boussinesq = true,
        include_turbulence = false,              # Laminar flow
        coupling_strategy = :segregated,         # Segregated coupling
        max_coupling_iterations = 8,
        coupling_tolerance = 1e-5,
        under_relaxation_factors = Dict(
            :momentum => 0.7,
            :pressure => 0.3,
            :temperature => 0.8,
            :turbulence => 0.5
        ),
        adaptive_relaxation = true,
        convergence_monitoring = verbose
    )
    
    if verbose
        println("Starting coupled simulation...")
        println("Expected convection time scale: t_conv ≈ $(L^2/α) s")
        println()
    end
    
    # Run coupled simulation
    result = solve_coupled_flow_heat_2d(ns_options, boussinesq_options, 
                                       heat_options, multiphysics_options)
    
    if verbose
        print_natural_convection_results(result, Ra, Pr)
    end
    
    return result
end

"""
    create_rayleigh_benard_bc(T_ref, ΔT) -> Dict{Symbol,Any}

Create boundary conditions for Rayleigh-Bénard convection.

# Arguments
- `T_ref::Float64`: Reference temperature (K)
- `ΔT::Float64`: Temperature difference between hot and cold walls (K)

# Returns
- `Dict{Symbol,Any}`: Boundary condition dictionary
"""
function create_rayleigh_benard_bc(T_ref::Float64, ΔT::Float64)
    T_hot = T_ref + ΔT/2    # Hot bottom wall
    T_cold = T_ref - ΔT/2   # Cold top wall
    
    return Dict{Symbol,Any}(
        :bottom => Dict(:type => :dirichlet, :value => T_hot),   # Heated bottom
        :top    => Dict(:type => :dirichlet, :value => T_cold),  # Cooled top
        :left   => Dict(:type => :neumann, :value => 0.0),       # Adiabatic left
        :right  => Dict(:type => :neumann, :value => 0.0)        # Adiabatic right
    )
end

"""
    print_natural_convection_results(result, Ra, Pr)

Print comprehensive results from natural convection simulation.
"""
function print_natural_convection_results(result::MultiPhysicsResult, Ra::Float64, Pr::Float64)
    println("\n" * "="^60)
    println("NATURAL CONVECTION SIMULATION RESULTS")
    println("="^60)
    
    # Simulation diagnostics
    println("Simulation Status:")
    println("  Converged: $(result.converged)")
    println("  Time steps: $(result.iterations)")
    println("  Total coupling iterations: $(sum(result.coupling_iterations))")
    println("  Average coupling iterations per step: $(mean(result.coupling_iterations))")
    println("  Solve time: $(result.solve_time) seconds")
    println("  Final residual: $(result.residual_norm)")
    println()
    
    # Physics diagnostics
    diag = result.physics_diagnostics
    println("Heat Transfer Results:")
    println("  Computed Rayleigh number: Ra = $(get(diag, \"rayleigh_number\", \"N/A\"))")
    println("  Average Nusselt number: Nu = $(get(diag, \"nusselt_number\", \"N/A\"))")
    println("  Average temperature: ⟨T⟩ = $(get(diag, \"average_temperature\", \"N/A\")) K")
    println("  Maximum velocity: |u|_max = $(get(diag, \"max_velocity\", \"N/A\")) m/s")
    println()
    
    # Flow field characteristics
    if !isnothing(result.u) && !isnothing(result.v)
        u_max = maximum(abs.(result.u))
        v_max = maximum(abs.(result.v))
        velocity_magnitude = sqrt.(result.u.^2 + result.v.^2)
        u_rms = sqrt(mean(result.u.^2))
        v_rms = sqrt(mean(result.v.^2))
        
        println("Flow Field Statistics:")
        println("  Max |u| velocity: $(u_max) m/s")
        println("  Max |v| velocity: $(v_max) m/s")
        println("  Max velocity magnitude: $(maximum(velocity_magnitude)) m/s")
        println("  RMS u velocity: $(u_rms) m/s")
        println("  RMS v velocity: $(v_rms) m/s")
        println()
    end
    
    # Temperature field characteristics
    if !isnothing(result.T)
        T_min = minimum(result.T)
        T_max = maximum(result.T)
        T_mean = mean(result.T)
        T_std = std(result.T)
        
        println("Temperature Field Statistics:")
        println("  Temperature range: [$(T_min), $(T_max)] K")
        println("  Mean temperature: $(T_mean) K")
        println("  Temperature std dev: $(T_std) K")
        println()
    end
    
    # Buoyancy effects
    if !isnothing(result.buoyancy_force)
        max_buoyancy = maximum(abs.(result.buoyancy_force))
        buoyancy_work = get(diag, "total_buoyancy_work", "N/A")
        
        println("Buoyancy Effects:")
        println("  Max buoyancy force: $(max_buoyancy) N/m³")
        println("  Total buoyancy work: $(buoyancy_work) W/m³")
        println()
    end
    
    # Dimensionless numbers
    println("Dimensionless Analysis:")
    println("  Input Rayleigh number: Ra = $(Ra)")
    println("  Input Prandtl number: Pr = $(Pr)")
    
    # Theoretical predictions
    if Ra < 1708
        println("  Regime: Conduction dominated (Ra < Ra_critical ≈ 1708)")
        println("  Expected Nu ≈ 1 (pure conduction)")
    elseif Ra < 1e5
        println("  Regime: Steady convection (1708 < Ra < 10⁵)")
        println("  Expected Nu ≈ 0.54 Ra^0.25 = $(0.54 * Ra^0.25)")
    else
        println("  Regime: Turbulent convection (Ra > 10⁵)")
        println("  Expected Nu ≈ 0.1 Ra^0.31 = $(0.1 * Ra^0.31)")
    end
    
    println("\n" * "="^60)
end

"""
    run_parametric_study(; Ra_range=[1e4, 5e4, 1e5, 5e5], Pr=0.71, resolution=24)

Run parametric study varying Rayleigh number.
"""
function run_parametric_study(; Ra_range::Vector{Float64}=[1e4, 5e4, 1e5, 5e5], 
                              Pr::Float64=0.71, resolution::Int=24)
    
    println("🔬 Natural Convection Parametric Study")
    println("=====================================")
    println("Rayleigh number range: $(Ra_range)")
    println("Prandtl number: Pr = $(Pr)")
    println("Grid resolution: $(resolution)×$(resolution)")
    println()
    
    results = []
    
    for (i, Ra) in enumerate(Ra_range)
        println("Case $(i)/$(length(Ra_range)): Ra = $(Ra)")
        println("-" * "="^40)
        
        result = run_natural_convection_demo(Ra=Ra, Pr=Pr, resolution=resolution, verbose=false)
        push!(results, result)
        
        # Print key results
        Nu = get(result.physics_diagnostics, "nusselt_number", NaN)
        max_vel = get(result.physics_diagnostics, "max_velocity", NaN)
        println("  Nusselt number: Nu = $(Nu)")
        println("  Max velocity: |u|_max = $(max_vel) m/s")
        println("  Converged: $(result.converged)")
        println()
    end
    
    # Summary comparison
    println("PARAMETRIC STUDY SUMMARY")
    println("="^50)
    println("Ra\\t\\tNu\\t\\tMax Vel\\t\\tConverged")
    println("-"^50)
    
    for (i, (Ra, result)) in enumerate(zip(Ra_range, results))
        Nu = get(result.physics_diagnostics, "nusselt_number", NaN)
        max_vel = get(result.physics_diagnostics, "max_velocity", NaN)
        converged = result.converged
        
        @printf "%.1e\\t%.3f\\t\\t%.3e\\t%s\\n" Ra Nu max_vel converged
    end
    
    return results
end

"""
    demonstrate_advanced_natural_convection()

Demonstrate advanced natural convection features including time-dependent heating.
"""
function demonstrate_advanced_natural_convection()
    println("🌡️ Advanced Natural Convection Demo")
    println("==================================")
    println("Features:")
    println("- Time-dependent boundary conditions")
    println("- Temperature-dependent properties")
    println("- Advanced thermal boundary conditions")
    println()
    
    # Create material properties with temperature dependence
    air_properties = create_air_properties(model=:temperature_dependent)
    
    # Base solver options
    ns_options = NSOptions(
        N = 4, nu = 1.5e-5, cfl = 0.25, tfinal = 10.0,
        dim = 2, n_block = 4, max_steps = 2000,
        tol = 1e-6, verbose = true
    )
    
    # Boussinesq options with realistic air properties
    boussinesq_options = BoussinesqOptions(
        reference_temperature = 300.0,
        thermal_expansion = 1/300.0,
        gravitational_acceleration = [0.0, -9.81],
        density_reference = 1.177,  # Air density at 300K
        prandtl_number = 0.71
    )
    
    # Advanced thermal boundary conditions
    boundary_system = create_thermal_boundary_system(thermal_conductivity=0.026)
    
    # Time-dependent heating on bottom wall
    time_dependent_heating = (x, y, t) -> 320.0 + 5.0 * sin(0.5 * t)
    hot_wall_bc = DirichletTemperature(
        temperature = time_dependent_heating,
        time_dependent = true
    )
    
    # Convective cooling on top wall
    convective_cooling_bc = ConvectiveBoundary(
        heat_transfer_coefficient = 25.0,  # W/m²·K
        ambient_temperature = 280.0,       # K
        time_dependent = false
    )
    
    # Add boundary conditions
    add_boundary_condition!(boundary_system, :bottom, hot_wall_bc)
    add_boundary_condition!(boundary_system, :top, convective_cooling_bc)
    
    # Heat transfer options with advanced boundary conditions
    heat_options = HeatTransferOptions(
        thermal_conductivity = 0.026,      # Air thermal conductivity
        thermal_diffusivity = 2.2e-5,     # Air thermal diffusivity
        heat_capacity = 1007.0,           # Air specific heat
        density = 1.177,                  # Air density
        boundary_conditions = Dict{Symbol,Any}(
            :bottom => Dict(:type => :dirichlet, :value => time_dependent_heating),
            :top => Dict(:type => :convective, 
                        :heat_transfer_coefficient => 25.0,
                        :ambient_temperature => 280.0),
            :left => Dict(:type => :neumann, :value => 0.0),
            :right => Dict(:type => :neumann, :value => 0.0)
        ),
        temperature_dependent_properties = true,
        time_integration = :implicit
    )
    
    # Multi-physics options
    multiphysics_options = MultiPhysicsOptions(
        coupling_strategy = :segregated,
        max_coupling_iterations = 10,
        coupling_tolerance = 1e-6
    )
    
    println("Running advanced natural convection simulation...")
    result = solve_coupled_flow_heat_2d(ns_options, boussinesq_options, 
                                       heat_options, multiphysics_options)
    
    print_natural_convection_results(result, 1e5, 0.71)
    
    return result
end

# Main execution
if abspath(PROGRAM_FILE) == @__FILE__
    println("Natural Convection Demo Examples")
    println("==============================")
    println()
    
    println("1. Basic Rayleigh-Bénard convection")
    basic_result = run_natural_convection_demo(Ra=1e5, verbose=true)
    
    println("\\n\\n2. Parametric study")
    parametric_results = run_parametric_study()
    
    println("\\n\\n3. Advanced features demonstration")
    advanced_result = demonstrate_advanced_natural_convection()
    
    println("\\n🎉 All natural convection demos completed successfully!")
end